O’Hara on GitHub


1 Summary

Examine biodiversity risk by taxonomic group.

2 Data sources

  • IUCN species API: IUCN. (2019). The IUCN Red List of Threatened Species. Version 2019-2.
  • IUCN species shapefiles: IUCN. (2019). The IUCN Red List of Threatened Species. Version 2019-2. Retrieved August 2019, from http://www.iucnredlist.org
  • BirdLife International shapefiles: BirdLife International and Handbook of the Birds of the World. (2018). Bird species distribution maps of the world. Version 7.0. Available at http://datazone.birdlife.org/species/requestdis
  • World Database on Protected Areas: IUCN, & UNEP-WCMC. (2018, June). The World Database on Protected Areas (WDPA). Retrieved June 9, 2018, from Cambridge, UK: UNEP-WCMC website: www.protectedplanet.net
  • Marine Ecoregions of the World: Spalding, M. D., Fox, H. E., Allen, G. R., Davidson, N., Ferdaña, Z. A., Finlayson, M. A. X., … others. (2007). Marine ecoregions of the world: A bioregionalization of coastal and shelf areas. BioScience, 57(7), 573–583.

3 Methods

3.1 set up taxonomic groups

3.1.1 maps of rr-risk by taxa

reload <- TRUE ### reload applies to building rasters, not the plot

rast_base <- raster::raster(file.path(dir_spatial, 'cell_id_mol.tif'))

land_poly <- sf::read_sf(file.path(dir_spatial, 'ne_10m_land', 
                                   'ne_10m_land_no_casp.shp')) %>%
  st_transform(crs(rast_base))

taxa_gps <- c("Cnidaria",
              "Tracheophyta",
              "Chordata: Actinopterygii",
              "Chordata: Mammalia",
              "Chordata: Chondrichthyes",
              "Chordata: Aves",
              "Chordata: Reptilia",
              "Echinodermata")

map_list <- vector('list', length = length(taxa_gps)) %>%
  setNames(taxa_gps)

for(i in seq_along(taxa_gps)) { ### i <- 6
    
  taxon_gp <- taxa_gps[i]
  taxon_txt <- taxon_gp %>% 
    tolower() %>% 
    str_replace_all('chordata: |[^a-z]+', '')

  taxa_info <- taxa_sums %>%
    filter(taxon == taxon_gp)
  
  taxa_rast_file <- file.path(dir_o_anx, 'taxa_rasts_mol_2019', 
                              paste0(tolower(taxon_txt), '_rr_comp.tif'))
  
  if(!file.exists(taxa_rast_file) | reload == TRUE) {
  
    cat_msg('Building: ', taxa_rast_file)

    rr_weighted_files <- taxa_info$rr_weighted %>% unique()
    
    taxon_df_rr <- parallel::mclapply(rr_weighted_files, mc.cores = 12,
        FUN = function (x) { ### x <- rr_weighted_files[1]
          read_csv(x, col_types = 'ddd___d___')
        }) %>% 
      bind_rows() %>%
      group_by(cell_id) %>%
      summarize(mean_risk_sum = 1/sum(sr_rr_risk) * sum(mean_risk * sr_rr_risk),
                pct_threat_sum = sum(sr_rr_threatened) / sum(sr_rr_risk)) %>%
      select(cell_id, mean_risk_sum, pct_threat_sum)
    
    mean_rr_rast <- raster::subs(x = rast_base, y = taxon_df_rr,
                                 by = 'cell_id', which = 'mean_risk_sum')
    
    raster::writeRaster(mean_rr_rast, filename = taxa_rast_file, overwrite = TRUE)
    
  } else {
    cat_msg('File exists: ', taxa_rast_file)
    mean_rr_rast <- raster(taxa_rast_file)
  }

  mean_rr_df <- mean_rr_rast %>% 
    aggregate(fact = 2) %>%
      ### should look OK for tiny maps?
    rasterToPoints() %>% 
    as.data.frame() %>%
    setNames(c('long', 'lat', 'value'))

  cat_msg('Building plot for ', taxon_gp)
  ### using geom_tile so I can add a size to pixels, thought it is
  ### significantly slower...
  x <- ggplot(mean_rr_df) +
    geom_tile(aes(long, lat, fill = value, color = value), 
              size = .5,
              show.legend = FALSE) +
    geom_sf(data = land_poly, aes(geometry = geometry), 
            fill = 'grey96', color = 'grey40', size = .10) +
    ggtheme_map(base_size = 9) +
    theme(plot.margin = unit(c(.05, .10, .05, .05), units = 'cm'),
          panel.background = element_rect(fill = 'grey80', color = NA),
          axis.title = element_text(face = 'bold', hjust = 1, 
                                    vjust = 0, size = 7,
                                    margin = margin(t = 0, r = 0, 
                                                    b = 0, l = .1, unit = 'cm')),
          axis.title.x = element_blank()) +
    coord_sf(datum = NA) + ### ditch graticules
    scale_fill_gradientn(colors = risk_cols, 
                         values = risk_vals, 
                         limits = c(0, 1),
                         labels = risk_lbls, 
                         breaks = risk_brks) +
    scale_color_gradientn(colors = risk_cols, 
                         values = risk_vals, 
                         limits = c(0, 1),
                         labels = risk_lbls, 
                         breaks = risk_brks) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    labs(y = taxon_gp)

  map_list[[i]] <- x
  
}